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O ■ ABSTRACT 

O' 
(N 

Qh' We describe the results of a simulation of coUisionless cold dark matter 

r^ . in a ACDM universe to examine the properties of objects collapsing at high 

ly-v . redshift {z = 10). We analyze the halos that form at these early times in 

this simulation and find that the results are similar to those of simulations of 
large scale structure formation at low redshift. In particular, we consider halo 
properties such as the mass function, density profile, halo shape, spin parameter. 



> 

lO . and angular momentum alignment with the minor axis. By understanding the 

properties of small scale structure formation at high redshift, we can better 



(N 



O ! understand the nature of the first structures in the universe, such as Population 

Q ! HI stars. 

O ' Subject headings: methods: n-body simulations — early universe 

^ 1. Introduction 

j^ ■ The formation of the first stars in the universe, also known as Population HI stars 

because of their lack of metals, is important for many reasons. These objects are responsible 
for the creation and dispersal of the first metals in the universe via Type II supernovae. The 
ultraviolet radiation from these first stars may also be partly responsible for the ionization 
of the intergalactic medium ([Haiman, Rees, fc Loeb 1997|) . These first objects are also the 
building blocks for the creation of larger structures such as galaxies in the bottom-up model 
of hierarchical structure formation. 

Various workers have studied the process of star formation in the early universe using 
numerical simulations of metal-free gas. Spherically symmetric simulations are useful 
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because they are computationally less intensive than three dimensional calculations and so 
can include more physics, such as detailed chemistry and cooling and even radiative transfer 
( [Haiman, 'I'houl fc Loeb 199(j| ; pmukai fc INishi 1998| ). However, to follow the collapse of 



matter to stellar densities, one needs a fully three dimensional approach. Several groups 
have carried out such simulations including the relevant chemistry ([Abel, Bryan, fc Norman 
2000t Promm, Coppi fc Larson 1999|) . 



These studies have focused on the evolution of individual density peaks, ignoring 
the effects of other matter in the universe. Tidal torques impart angular momentum 
to collapsing objects, affecting their subsequent evolution. Virialized halos can merge, 
changing the structure of objects. In this paper, we model a relatively large volume, 
1 h^^Mpc (comoving), that contains a substantial number of halos, in order to understand 
the nature of the first objects that collapse on a more statistical level. 

Since we are interested in the gravitational behavior of the matter, we can model the 
evolution with an N-body simulation. This is a valid approximation for a universe where 
collisionless cold dark matter dominates. In such a universe, the gas is coupled to the dark 
matter, so that gas falls into the potential wells of dark matter halos which may in turn 
lead to star formation. Thus, we can gain some understanding of the sites of the formation 
of the first stars via simulations of dark matter. 

Thus far, numerical simulations with dark matter have focused on the problem of large 
scale structure formation (e.g. pUfstathiou, et al. 1988| ; [Katz, Hernquist, fc Weinberg 1999| ; 



Kauffman, et al. 1999| ; |Lacey fc Cole 199^) . These simulations address the question of 



structure formation on the scale of galaxy clusters in order to understand the processes of 
galaxy formation and clustering. Here, we apply the analytic tools that these groups have 
developed to address the question of structure formation on a smaller scale. 

This paper describes the properties of collapsed objects in a numerical simulation at 
high redshift and small scales in a ACDM universe. These objects are some of the first 
non-linear structures in the universe. The non-linear evolution of the power spectrum 
has been addressed separately ( |Jang-Condell fc Hernquist 2000| ). In §2, we summarize 



the computational method used for the simulation; in §3, we describe the results of the 
simulation; and in §4 we summarize our results and discuss them in comparison to previous 
work with N-body simulations. 
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2. Method 

In this simulation, we model the evolution of dark matter in a periodic cube of 
size 1 h~^Mpc per edge in comoving coordinates. We adopt a ACDM cosmology, with 
Qm = 0.35, Q\ = 0.65, and h = 0.65, where the Hubble constant is Hq = lOOh km/s/Mpc. 
The power spectrum is normalized to as = 0.9, where as is the rms density variation 
smoothed with a top-hat filter of radius 8 h~^Mpc. For the cosmological model we consider, 
this normalization is roughly consistent with both the local abundance of rich clusters 
( [White, Efstathiou, fc Frenk 199.^ ) and fluctuations in the cosmic microwave background as 
observed by COBE ( |Bennett, et al. 1996| ). 

Initial conditions were generated using the analytic fit to the CDM power spectrum 
derived by Efstathiou, Bond, & White (1992) 

Bk 



P{k) = \5^ 



|2 



{1 + [ak + (6A;)3/2 + (cA;)2]^}2/^ ' ' ' 

where a = 6.4/r h'^Mpc, b = 3.0/T /i^^Mpc, c = 1.7/r /i^^Mpc, z/ = 1.13, and 5 is a 
normalization constant determined by erg. The shape parameter F is set to F = Qmh for 
a ACDM cosmology. This power spectrum is extrapolated to a starting redshift of 100 
assuming linear theory, and is then converted into spatial density fluctuations by first 
assigning random phases to the 6k and then taking the Fourier transform. The spatial 
density fluctuations are converted into particle positions and velocities using the Zel'dovich 
approximation. 

The code used for the simulation is PTreeSPH, a gravity treecode with smoothed particle 
hydrodynamics (SPH) designed to run on a parallel supercomputer (Dave, Dubinski, & 
Hernquist 1997). The SPH part of the code was unused in this simulation since gas was 
not included. The code uses a Barnes-Hut ([Barnes fc Hut 1986] , Pernquist 1987] ) algorithm 



for computational efficiency in calculating gravitational forces and a spline kernel for 
gravitational softening ( [Hernquist fc Katz 1989| ). The softening length chosen for this 
simulation was l/20th of the mean interparticle separation, or l/2560th of the box size. The 
simulation box contains 128^ particles, implying that each particle represents 4.6 x IO^Mq 
of dark matter. The simulation was run on a four-processor Beowulf-type cluster of PCs at 
the Harvard-Smithsonian Center for Astrophysics. 



3. Results 

The particle distribution at the final output of redshift z = 10 is shown in Fig. |1[ The 
particles are displayed in projection along the x-axis and the colors indicate the particle 
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densities. Already at this high redshift there is evidence of structure formation in the form 
of clumps and filaments of high particle density. The structures are similar to those seen in 
large scale simulations evolved to much lower redshift. 

As we can see in Fig. ^ the centers of dense knots have densities of ^ lOOOp. Density 
evolves with redshift as a~^ = (1 + z)^, so at 2; = 10, lOOOp = 1000a~^Qi,pcrit where Qf, is 
the baryon density of the universe in units of the critical density, and pcrit is the present day 
critical density. Constraints on light element production during big bang nucleosynthesis 
indicate that 0.01 ^ Qbh'^ ^ 0.015 (e.g. Peacock 1999| , §9.5). Thus, 



lOOOp ^ 3 X 10"^^ g cm"^ (2) 

which is ~ 0.2 cm~^ in hydrogen atoms. This is comparable to the density of neutral 
hydrogen in the present day ISM, which is ~ 1 cm~^ (e.g. [Spitzer 1978|, §1.1). We can see 
that the densest regions in the simulation will be similar in density to the present-day ISM, 
and so are potential sites for star formation. 



3.1. Finding halos 

The particle distribution was analyzed using a program called skid to determine the 
positions, masses, and sizes of collapsed halos. A description of skid can be found at 
tittp : //www-hpcc . astro . Washington . edu/tools/SKID/| . 



The basic algorithm that skid uses is as follows: 

1. Calculate densities, and consider only those above a user specified density threshold. 
These are called the moving particles. 

2. Slide the moving particles along the density gradient toward higher density. 

3. Continue moving particles until all the particles stop moving and are localized in high 
density regions of a user specified size (eps). 

4. Group together these localized particles using the friends-of-friends method with a 
linking-length of eps. 

5. Reject groups with less than a user specified minimum number of particles. 

6. Particles which are not bound to their group are removed from the group. 

7. Reject groups with less than minimum number of particles. 



-5- 




Fig. 1. — Simulation output at redshift z = 10. The particles are shown in 
projection along the x-axis. Color indicates log density in units of p. The size 
of the box is 1 h~^Mpc (comoving) per side. (This figure is also available at 



tittp : //cf a-www . harvard . edu/~h-j ang/research/prettypic . gif 
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We used a density threshold of 200p, a hnking-length (eps) of l/640th of the box 
size or 1/5 the mean interparticle spacing, and a minimum halo size of eight particles. 
Using the simulation output at a redshift of ^ = 10, this resulted in 2881 halos, with the 
largest halo consisting of 8475 particles, corresponding to a mass of 4 x IO^Mq and the 
smallest with eight particles, corresponding to 4 x IO^Mq. The density cutoff was chosen 
somewhat arbitrarily, and changing its value does not significantly affect the overall results 
presented in this paper. For example, changing the density threshold from 200p to 86p did 
not alter the ensemble properties of the halos. Some of the halos identified by skid were 
subsequently rejected from the sample for being unbound, as described in § p.3| . We have 
analyzed various properties of these halos and describe our results below. 



3.2. Mass function 



The distribution of halos masses can be expressed in terms of the mass function, /(M), 
where f{M)dM is the number density of halos of with mass between M and M + dM. 
An analytic prediction for the mass function can be obtained using the Press-Schechter 
formalism ( [Press fc Schechter 1973| ). 



The Press-Schechter formalism states that the fraction of the universe condensed into 
objects of mass > M is 



F(> M) = 1 - erf 



(3) 



,V2a{M)J 

where a{M) is the rms density fluctuation smoothed over spheres of mass M, and 5c is the 
critical overdensity. The critical overdensity is defined as follows: an object of density p in 
a universe of average density p is collapsed when 



s = f^^>s. 



(4) 



We take 6c = 1.69, the canonical value. The mass function then becomes 



/(M) dM 



dF 



dM 



U^ = \\--p 

M V vr a 



dXwa 



dlnM 



exp 



5^. \ dM 



2^2/ M2 



(5) 



Figure ^ shows the actual mass function of halos in the simulation compared to the 
predictions of the Press-Schecter formalism. Halos with positive binding energy were 



omitted from the calculation of the mass function as explained in Section ^^. The results 
are in remarkably close agreement with the Press-Schechter prediction. 
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Fig. 2. — Mass function of halos. The points represent the simulation data and the sohd hne 
is the Press-Schechter prediction, with 5c = 1.69. 



3.3. Halo energies 

When skid calculates lialos from the particle distributions, it removes unbound 
particles from the halos, so in principle, each halo should be gravitationally bound. The 
binding energy of particle i to the rest of the halo is 



Ei = rrii 



i + E^A^ 



(6) 



where m^ and rnj are the masses of particle i and particle j respectively, Vi is the velocity 
of particle i with respect to the halo center of mass velocity, and $j(rjj) is the gravitational 
potential between particles i and j. The gravitational potential has the general form 
$j(rjj) = —Gmjf{rij), which admits a softened potential for particles pairs with small rjj 
( IHernquist fc Katz 1989]) . 

The requirement for binding is that E^ < 0. The binding energy for the halo as a whole 
is the sum of its gravitational potential energy and its kinetic energy. Thus, 

E = K + W = ^J2 m^vl - ^ E E Grn,rn,f{n,). (7) 

i i i^j 

The factor of | in the potential energy is to account for double- counting the interactions 
between pairs of particles. Equation (^ can be rewritten as 

^ = E ^^ + ^ E E Gm,m,fiu,) = Y.E^ + \W\ (8) 

i i i^j i 

which is greater than the sum of the individual binding energies of the particles. In other 
words, it is possible to have a halo where each particle is bound to the sum of all the other 
particles, but the halo as a whole is unbound. So, after skid has calculated the halos, the 
total binding energy of each halo is checked, and discarded from the sample if it is unbound. 
Out of the total of 2881 halos, 380 were rejected in this way, leaving 2501 bound halos. 

The distribution in size of the 380 unbound halos is plotted in Figure 0. Note that the 
unbound halos are all quite small, the largest containing only 25 particles. This shows that 
our method of using skid to calculate halos is fairly robust, failing only at small halo sizes, 
which are intrinsically unstable to small perturbations. In the analysis of the remaining 
halo properties, calculations are done either as a function of mass, or neglecting all halos 
with less than 50 particles. Thus, the removal of these halos from the sample will not 
significantly affect our results. 

Figure ^ is a plot of the energies of the halos versus their masses. The absolute values 
of the potential energy are plotted as squares and the absolute values of the total binding 
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Fig. 3. — Distribution of unbound halos with respect to halo size in terms of particle number. 
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Fig. 4. — Energies of the halos versus mass. Potential energy is plotted as squares, total 
binding energy is plotted as crosses. The solid line is proportional to M^/^. 
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energy are plotted as crosses. The energies all seem to follow a power law oi E oc M^/^, as 
indicated by the solid line. This relation can be obtained analytically by supposing that 
the halos are the result of the collapse of spheres of uniform density. This is a reasonable 
assumption since the universe was initially very nearly uniform in density. If the density is 
constant, then the mass of a sphere of radius r is 

M{r) = ^r\ (9) 

The equation for the gravitational potential energy energy of a uniform sphere of density p 
and mass M is 



f 
W = (^{r)pdPr 

Jo 

= / —^pA-nrdr 

Jo r 



2 



3GM 
5R 

iGM=/3fl!V)'". (10) 



Assuming that the halo is virialized, we have 

E = K + W = -^W + W = ^W. (11) 

Thus, the total energy should obey the same scaling. 

The dotted line represents 1/2 the energy of the solid line. The total energies of the 
halos are all below this line, indicating that E > ^W in general. This means that the 
halos have not yet become virialized. This is also shown in Fig. ^, which shows the ratio 
between the kinetic energy and the potential energy. The solid line marks K = ^\W\, which 
would indicate a virialized halo. Nearly all the halos have too high a kinetic energy for 
virialization. 



3.4. Gas temperature 

Assuming that the matter in the universe is dominated by cold dark matter, then the 
gas should trace the distribution of the dark matter. In particular, gas should fall into the 
potential wells of the dark matter halos and become shock heated to the virial temperature 
of the halo. We can estimate this temperature by assuming that the mean-square velocity 
of the gas is the same as the dark matter. 
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Fig. 5. — Ratio of kinetic to potential energy. The solid line indicates K 
holds for virialized objects. 
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For an ideal gas, the temperature is related to the mean-square velocity by 

^' = , 12 

/imp 

where h is Boltzmann's constant, T is the temperature in Kelvin, m^ is the proton mass, 
and \i is the mean molecular weight. This temperature, as a function of /x, is plotted versus 
mass for each halo in Fig. p. The value of /x will depend on the composition and ionization 
state of the gas. 



As discussed in section p.3| , the kinetic energy of the halo follows a power law of 
K oc M^/'^. Since K = ^M{v'^), the temperature should follow the relation 

T ~ M^/\ (13) 

This is indicated by the solid line in Fig. |[ The temperatures follow this power law relation 
fairly well. 



3.5. Halo profiles 

We calculated the density profiles of the thirty most massive halos. The spherically 
averaged density as a function of radius for these halos, with the potential minimum as the 
center, is plotted in Figs. |^ and p. The potential minimum was chosen over the center of 
mass because the potential minimum is more likely to be at the point of highest density. 

Possible analytic fits to dark matter halo profiles are those due to Hernquist (1990) 
and Navarro, Frenk & White (1997, henceforth NFW). Avila- Reese, et al. (1999) proposed 
the following general form for the fitting formulae: 



P(0 = . '' ,,-1 (14) 



Ji(l + ^ 

where po and r^ are scalings for the density and radius, respectively, and /3 is a power-law 
index parameter explained as follows. When /3 = 4 this corresponds to a Hernquist profile, 
and when when (3 = 3 this corresponds to an NFW profile. We can also allow j3 to he a. free 
parameter. The asymptotic behavior of the profiles is 

f r-\ r/rs < 1 . . 

P ^ [ r-^ r/r, » 1. ^'^^ 

Thus, (3 parametrizes the profile shape at at radii much larger than the scale length. 
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Fig. 6. — Virial temperature versus mass. The solid line is proportional to M^/^. 
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The fits to the NFW and Hernquist profiles are indicated in Figs. |^ and |^ by dotted 
and short dashed fines, respectively. A third model, allowing /? to be a free parameter in 
the fit, is indicated by a long dashed line. This fitted (3 is displayed at the bottom of each 
profile plot. The scale radius r^ for each of these three profiles is indicated by filled, fat, 
and thin triangles, respectively. 

All three profiles are good fits to the data, and are indistiguishable. This is because 
the power- law slopes for all the halos are fairly close to —2. Since each profile is a smoothly 
varying function, the power-law slope also varies smoothly from —1 to —j3. One can easily 
fit the part of the profile with power-law slope —2 to each halo profile. Note also that 
(3 ~ 3.5 in all cases, intermediate to the NFW and Hernquist profile expectations. The close 
agreement of the fits can also be explained by the relatively large values of r^ with respect 
to the sizes of the halos. The halos are only a few scale radii in size, so we are well below 
the regime where the halo profiles behave as p oc r~^. 



3.6. Spin parameter 

The spin parameter of a bound object is defined as 

where J is the angular momentum, E = K + W is the binding energy, G is the gravitational 
constant, and M is the mass. This is a dimensionless parameter that quantifies the amount 
of rotational support an object has. A value of A ~ 1 corresponds to nearly full rotational 
support, and is typical of spiral galaxies, while A ~ 0.05 is typical of elliptical galaxies, 
which are supported by velocity dispersion. Figure ^ is a plot of spin parameter versus mass 
for the halos in the simulation. 

Mo, Mao & White (1998) found that the distribution of spin parameters of dark matter 
halos in galaxy formation simulations approximates a log-normal distribution, that is 



p(A) d\ = — exp 
v27rcr 



ln^(A/A) 
2a2 



T- '-) 



with fitting parameters of A = 0.05 and ax = 0.5. The log-normal distribution is also a 
good fit for the data presented in this paper, as shown in Fig. ITO. The simulation data is 



displayed as squares, and the log-normal fit to it is the solid line. The fitting parameters are 
A = 0.043 and cxx = 0.53, which is consistent with the results of Mo, Mao & White (1998). 

The smallest halos in the simulation contain only eight particles. For a halo with so 
few particles, the spin parameter is not a particularly meaningful quantity. Since the most 
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Fig. 7. — Density profiles of halos. Radius is in physical parsecs. Dotted line is NFW profile 
fit, short dashed line is Hernquist profile fit, and long dashed line is the fit allowing /5 to be a 
free parameter. Triangles mark the r^ for each fitted profile - solid triangle is NFW, hollow 
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Fig. 9. — Spin parameter versus mass. 
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numerous halos in the simulation are the smallest ones, the distribution of A calculated 
above may not be particularly meaningful either. We recalculated the distribution of A for 
halos containing at least 50 particles to determine the effect of excluding small halos. This 
is plotted as diamonds in Fig. lU, and the log- normal fit to it is the dashed line. The fitting 
parameters for these 389 halos is A = 0.033 and ax = 0.52, indicating that the larger halos 
have systematically smaller spin parameters than smaller halos. 



3.7. Halo shapes 

A typical dark matter halo is not actually spherical, but is triaxial due to anisotropy 
in its velocity dispersion. Thus, the shape of a halo can be quantified by finding its best 
fitting ellipsoid. The three axes of the best fitting ellipsoid will be referred to as the major, 
intermediate, and minor axes, in decreasing order of size. 



3.7.1. Method of calculating ellipsoids 

The best-fitting ellipsoid can be found by using the moment of inertia tensor, which is 
defined as 

E rriiiyf + 4) - E miXiVi - E rriiXiZi 

i i i 
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i i i 
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where r^ = Jxl -|- yf + zf is the distance of the zth particle to the center of the distribution 
of particles. The moment of inertia I about an axis through the center of the distribution 
in the direction of the unit vector n is given by 



n ■ I ■ n. 



(19) 



Another property of the moment of inertia tensor is that the angular momentum L can be 
expressed as 

L = I ■ cj. (20) 
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Fig. 10. — Distribution of spins. The squares and diamonds show the normahzed distribution 
for all the halos and halos with at least 50 particles, respectively. The solid line and the 
dashed line show the log-normal fit to the respective data sets. 
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Thus, the eigenvectors correspond to the axes about which the angular velocity and angular 
momentum are aligned and the eigenvalues are the corresponding moments of inertia. 

Consider an ellipsoid of constant density centered at the origin with its axes along the 
X-, y-, and 2;- axes defined by 

(21) 



x^ u^ z"^ 
h — H < 1. 



a" 



62 



c" 



By symmetry, the moment of inertia tensor is already diagonalized, and the diagonal 
elements are the eigenvalues. It can be shown (through some rather tedious integration) 
that the diagonal elements are 
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(22) 



where M is the total mass of the ellipsoid. Using equation (|T8|), we now have 
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which reduces to 
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Thus, the axes of the best fitting ellipsoid can found by calculating the tensor 
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and finding its eigenvalues and eigenvectors. 

The simplest way of determining the halo shape is to use the center of mass as the 
center of the distribution and to sum over all the particles in the halo to calculate the tensor 
Map. However, some halos found by the skid program include satellite halos or consist of 
two or more groups connected by a thin bridge. The center of mass for these halos may not 
actually lie in the center of the halo. In addition, ellipsoids calculated about the center of 
mass have systematically higher axial ratios - that is, they are more spheroidal in shape. 
This can be explained by the fact that skid tends to calculate halos within a spherical 
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volume, regardless of the intrinsic shape of the halo. Thus, the calculated shape of the halo 
is more rounded out, so to speak. 

A better approximation is to use the potential minimum as the center of the distribution 
and to sum over the inner part of the halo. We used an iterative method to fit an ellipsoid 
to half the mass of the halo. The procedure begins by finding the sphere centered at the 
potential minimum that contains half the mass of the halo. These particles are used to 
calculate an initial guess for the ellipsoid axes. Keeping the orientation and axial ratios of 
the ellipsoid fixed, a new ellipsoid is calculated which contains half the mass of the halo. 
The particles within this ellipsoid are used to calculate a new Ma/s and a new guess for the 
ellipsoid is calculated. This procedure is repeated until the values of the axes converge. 
Dubinski and Carlberg (1991) employed a similar method to calculate halo shapes, but 
using particles within a fixed distance of the halo center rather than a fixed fraction of 
the mass. In order to leave out halos that are too small to have their shapes accurately 
calculated, we chose a lower cutoff to the halo size of 50 particles. There were 389 halos of 
this size. 



3.7.2. Ellipticities and triaxiality 

Using the iterative procedure described above, we calculated the magnitude and 
orientation of the principal axes of each halo. Henceforth, we shall refer to the lengths of 
the major, intermediate, and minor axes as a, b, and c, respectively. The ellipticities of the 
halos, 

ei = 1 — b/a, 
£2 = 1- c/b, 

are useful measures of the shapes. A perfectly spherical halo would have a = b = c, hence 
ei = €2 = 0. An oblate halo would have the larger two axes equal to each other (a = 6), 
hence ei = and e2 > 0. A prolate halo would have the smaller two axes equal to each 
other (6 = c), hence ei > and €2 = 0. In general, however, a halo will be triaxial, meaning 
that there are no equalities between axis lengths. 

Figure |ll| is a plot showing the ellipticities of the halos. The diagonal line shows the 
division between prolate and oblate halos - prolate halos lie below the line, and oblate 
halos lie above it. By this criterion, there are 212 prolate halos and 177 oblate halos. 
Warren et al. (1992) also found more prolate than oblate halos in their simulations of dark 
matter halos at zero redshift at galaxy size scales (the smallest mass they considered was 
3 X IO^Mq). They found this result to be independent of the initial power spectrum used, 
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so it is reassuring that we obtain the same results for our simulation, although we use a 
different power spectrum and mass range. 

We can also compare the prolateness/oblateness of halos by use of the triaxiality 
parameter T, which is defined as 

T = -. 26 

A halo that is purely prolate has 6 = c, so T = 1. A halo that is purely oblate has a = 6, so 



T = 0. The dotted lines in Fig. 11 represent contours of constant T. 



Figure |T2| is a plot of T versus halo mass. There does not appear to be any strong 
correlation between mass and triaxiality. However, there is a tendency at all masses for the 
triaxiality to be close to 1. Figure |13| shows the distribution of the triaxiality parameter 
for all halos with at least 50 particles. The halos tend toward high values of T, another 
indication that prolate halos dominate. These two figures are qualitatively similar to 
Figs. 9a and 8 in Warren et al. (1992), indicating that our results are similar to theirs, 
despite the difference in mass scales and redshift. 



3.7.3. Angular momentum misalignment 

Now that we have calculated the axes of the best-fit ellipsoids of the halos, we can 
consider the orientation with respect to the angular momentum. Since the ellipsoids were 
fit using half the mass of the halo, we also recalculate the angular momenta using the 
positions and velocities of the same particles that were used to calculate the ellipsoids. We 
find the cosine of the angles between the angular momentum and each of the principal axes 
by taking dot products. 

Figure |l^ shows the distribution of these angle cosines. The dotted, dashed, and solid 
lines show the distribution of the angle cosines between the major, intermediate, and minor 



axes of the halo, respectively. As we can see in Fig. 0, the angular momentum tends to 
be aligned with the minor axis of halos. This is consistent with observations of elliptical 
galaxies which indicate that the angular momentum vector tends to align with the the 
minor axis ( |Franx, et al. 199TD . 



Kinematically, the angular momentum should align with either the major or minor axis 
because particles in a triaxial potential admit tube orbits only about these axes. In other 
words, particles may only circulate about either the major or the minor axis, causing the 
angular momentum of the system as a whole to align with these axes ( pinney fc Tremaine| 



1994| ). However, most halos do not show alignment with the major axis. In fact, many have 
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Fig. 11. — Ellipticities of halos of size > 50 particles, defined as ei = 1 — b/a and £2 = 1 — c/b. 
The dashed line shows the division between prolate and oblate halos, and the dotted lines 
show contours of constant triaxiality. 
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Fig. 12. — Triaxiality of halos with > 50 particles versus mass. 
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Fig. 13. — Distribution of halo triaxiality parameter T, defined as T = (a^ — &^)/(a^ — c^ 
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Fig. 14. — Distribution of angle cosines between angular momentum vector and ellipsoid 
axes for halos with > 50 particles. The dotted (dashed, solid) line shows the cosine of the 
angle between the major (intermediate, minor) axis of the halo. 
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their angular moinentuiii perpendicular to their major axis. 

In addition, there are a number of halos whose angular momentum is aligned with 
the intermediate axis, even though this configuration is kinematically unstable. Warren et 
al. (1992) also observe this in their dark matter simulations, and suggest that this is caused 
by a long time scale for the realignment of orbits. 



4. Summary 

The major conclusion of that paper is that this N-body simulation of small scale 
structure formation at high redshift is similar to N-body simulations of large scale structure 
formation at at low redshift. The details can be summarized as follows: 

1. The mass function of the halos can be well approximated by the Press-Schechter 
formalism. 

2. The profile of the halos can be modeled equally well by either the NFW and Hernquist 
profiles. 

3. The shapes of the halos are generally triaxial, with a tendency toward prolateness. 

4. The average spin parameter of the halos is about 0.04. 

5. The angular momentum tends to align with the minor axis most often, and favors 
alignment with the intermediate axis over the major axis. 

In fact, it would be surprising if the results were completely different from other N-body 
simulations, since they all model the gravitational collapse of objects from some primordial 
power spectrum of density fluctuations. The main difference is in the shape of the power 
spectrum, which changes as we go to smaller scales. This produces a corresponding change 
in the mass spectrum of objects that are collapsing. 

The Press-Schechter formalism predicts how the mass spectrum of collapsed objects 
depends on the initial power spectrum. We flnd that the mass function matches the 
predictions of Press-Schechter remarkably well, even at the low mass end. This is consistent 
with previous work with N-body simulations at low redshift and large scales which have 
also shown good agreement with the Press-Schechter formalism ( [Efstathiou, et al. 1988| ; 



Kauffman fc White 1993| ; |Lacey fc Cole 199^ ; [Valageas, Lacey fc Schaeffer 2000| ). 



The remaining halo properties described in this paper, including density proflle, halo 
shape, and spin parameter, all are consistent with previous work with N-body simulations 
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regardless of the power spectrum used. Some authors use scale invariant power spectra, 
but using a range of spectral indices: choices of —2, —1, and are typical ( |Cole fc Lacey 



1996| ; [Warren, et al. 19"^ ) . Some use a CDM spectrum or some variant ( [Mo, Mao, fc White 



1998 ). Note that at the high redshifts studied in this paper, the dynamical evolution of a 
ACDM universe is the same as that of a pure CDM universe, since the matter density varies 
as (1 + z)^ while the vacuum energy density remains constant. However, the spectral index 
approaches ~ —3 on small scales, so we effectively choose a different power spectrum by 
considering small scales. Although the initial power spectrum used in this paper is different 
from previous work, the overall results on halo properties remains the same. 

The shapes of the density profiles of the halos are consistent with both NFW and 
Hernquist profiles, which are based on studies of collapsed halos in cold dark matter 
simulations. Indeed, when the shape of the outer profile (3 is allowed to be a free parameter, 
we find that 3 < /? < 4, intermediate to the NFW and Hernquist profiles. 

The results for the spin parameter and shapes of halos are also consistent with 
large-scale structure simulations. The median and distribution of the spin parameter are 
both similar to the those found for large-scale simulations, i.e. A = 0.043, with a log-normal 
distribution ( |Mo, Mao, fc White 199^ ; Cole fc Lacey 19^ ). The predominance of prolate 



halo shapes also agrees with these other simulations ( [Warren, et al. 199^ ; [Cole fc Lacey[ 



1996[) , as well as the alignment of the angular momentum with the minor axis (jWeil & 



Hernquist 1994[ ; [Warren, et al. 1992[) . We also find that the angular momentum favors 



alignment with the intermediate axes over the major axis, as Warren et al. (1992) do. 

We can infer from this simulation that dark matter halos on small scales at high redshift 
behave very similarly to halos on large scales at low redshift. This can help us understand 
high-redshift star formation by providing information on the dark matter environments in 
which the first stars formed. Further study using simulations such as this may shed light on 
the IMF of the first stars, and their subsequent fate. By adding additional physics to the 
simulation, such as gas physics and radiative transfer, we can study the nature of the first 
objects that formed in the universe. 
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